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NOTATION 


E Expectation  (mean  value)  operator 

f^(w)  Wave  spectral  density  function 

R(t)  Covariance  function 

Re  Real  part 

S (w)  Wave  energy  spectrum 

T Wave  period 

t Time 

X(t)  Random  function  of  time 

Random  phase  angle  (discrete) 

e( to)  Random  phase  angle  (continuous) 

n(t)  Wave  surface  elevation 

2 

a Variance  of  wave  record 

T Time  lag 

<px(p)  Characteristic  function  of  X 

u)  Radian  frequency 

u>k  Discrete  frequency  component 
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ABSTRACT 


This  report  presents  a convenient  method  for  simu- 
lating time  histories  of  ocean  wave  surfaces  on  a digital 
computer.  The  simulation  model  employed  is  basically  the 
spectral  energy  wave  surface  representation  of  St.  Denis 
and  Pierson.  A random  selection  of  frequency  components 
within  each  frequency  interval  is  employed  to  increase 
the  wave  period  of  the  digital  process.  Analytical 
properties  of  the  resulting  wave  system  are  derived  and 
comparison  with  computed  wave  records  is  shown.  Detailed 
computer  program  documentation  is  also  provided. 


ADMINISTRATIVE  INFORMATION 

This  work  was  performed  under  the  General  Hydromechanics  Research  Pro- 
gram of  the  David  W.  Taylor  Naval  Ship  Research  and  Development  Center  and 
was  funded  by  the  Naval  Sea  Systems  Command  under  Program  Element  R02301, 
Project  61153N,  Task  Area  SR  0230101,  and  Work  Unit  1524-470. 

CONVERSION  TABLE 

Measurements  are  given  in  customary  units.  Their  conversion  to  metric 
equivalents  can  be  computed  using  the  following  values: 

1 ft/sec^  = 0.3048  m/s^ 

INTRODUCTION 

Time  histories  of  wave  surface  elevation  are  employed  extensively  in 
the  study  of  ship  performance  in  a seaway.  Such  time  series  can  be  used 
both  to  study  properties  of  the  wave  system  itself  and  to  provide  ex- 
citation for  the  various  ship  motions  and  propulsive  responses.  Current 
interest  in  the  time  domain  simulation  of  ship-wave  dynamic  systems  has  led 
to  the  need  for  digital  synthesis  of  wave  time  histories.  This  report 
presents  or.e  method  for  generating  such  time  histories  on  a digital 

computer.  The  representation  of  wave  surface  elevation  is  based  on  the 

1* 

random  phase  wave  energy  integral  model  used  by  St.  Denis  and  Pierson, 
with  two  modifications:  first,  the  wave  representation  is  interpreted  as 
a stochastic  integral  which  is  closely  related  to  the  Riemann-Stielt jes 
integral;  second,  the  frequency  values  for  determining  the  spectral 


*A  complete  listing  of  references  is  given  on  page  65. 
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amplitudes  and  phases  are  chosen  randomly  within  each  interval  in  order  to 
increase  indefinitely  the  period  of  the  resulting  digital  process. 


BACKGROUND 

It  is  commonly  accepted  in  linear  wave  theory  that  points  of  wave 

surface  elevation  follow  a normal  (Gaussian)  distribution.  This  Gaussian 

property  of  ocean  waves  was  demonstrated  in  a statistical  analysis  of 

2 

field  data  over  two  decades  ago  by  Rudnick  and  it  has  since  been  supported 
by  various  studies.  Various  past  investigations  have  dealt  with  the  sea 
surface  as  a stationary  Gaussian  process,  which  represents  approximately  a 
relatively  stable  sea  condition. 

Many  representations  of  stationary  Gaussian  processes  are  possible. 

For  example,  the  earliest  representation  was  through  the  use  of  Fourier 

series  expansions  having  normally  distributed  coefficients.  Such  models 

were  used  as  early  as  1910  by  Einstein  in  his  investigation  of  black-body 

3 

radiation  and  they  are  discussed  in  detail  by  Rice.  Fourier  series  models 
are  appealing  for  application  to  ocean  waves  because  of  their  nice 
physical  interpretation.  However,  they  have  the  unfortunate  disadvantage 
of  being  periodic.  Their  periods  can  be  increased  by  lowering  the  funda- 
mental frequency  of  oscillation  which  may  result  in  an  extremely  large 
number  of  frequency  components  in  order  to  cover  the  full  frequency  range 

4 

of  interest.  More  recently,  Levy  has  given  a different  representation, 
the  so  called  "canonical  representation"  of  Gaussian  processes,  in  terms 
of  integrated  Gaussian  white  noise.  This  model,  which  requires  independent 
Gaussian  random  variables  for  its  generation,  would  seem  to  involve  ex- 
tensive computational  time  in  practical  application.  Fourier  series 
expansions  containing  deterministic  coefficients  but  uniformly  distributed 
phase  angles  have  also  been  employed  to  obtain  periodic  processes  that  are 
asymptotically  normal  as  the  number  of  frequency  components  approaches 
infinity  (Grenander-Rosenblatt"*) . 

One  representation  of  a stationary  Gaussian  wave  surface  that  appears 
well-suited  for  practical  application  is  the  continuous  frequency  random 
phase  model  which  was  introduced  by  Levy**  and  which  was  proposed  as  a 
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model  for  ocean  waves  by  St.  Denis  and  Pierson.  This  model  yields  a 
nonperiodic  random  phase  representation  of  the  wave  surface.  The  St. 
Denis-Pierson  model  has  been  used  extensively  in  frequency  domain  in- 
vestigations of  ship  performance;  however,  it  appears  that  the  model  has 
not  been  explicitly  exploited  for  time  domain  prediction  or  simulation. 

This  report  implements  the  model  in  a time  domain  digital-computer 
simulation. 

THEORY 

The  random  phase  ocean  wave  model  may  be  represented  by  the  stochastic 
integral 


OO 


In  this  representation  n(t)  is  the  unidirectional  wave  surface  elevation. 
The  random  phase  angles,  e(io),  are  a function  of  frequency  and  are 
assumed  to  belong  to  an  uniform  distribution  over  the  interval  [0,2u). 

S(iu)  is  the  wave  energy  spectrum,  which  gives  the  mean  squared  value  of 
amplitude  associated  with  each  non-negative  frequency.  The  area  under  the 
spectrum  equals  the  variance  of  the  wave  record.  Thus,  S(w)  is  intimately 
related  to  the  sea  state.  It  is  assumed  that  the  sea  state  is  stationary 
so  that  the  energy  spectrum  will  be  independent  of  time. 

For  simplicity,  the  wave  surface  model  is  here  assumed  to  be  uni- 
directional and  spatially  fixed  at  the  origin  of  the  wave  field.  That  is, 
the  model  represents  surface  elevations  at  a point  in  space  and  considers 
only  long-crested  waves.  Extension  to  multidirectional  and  spatially 
dependent  models  is  straightforward. 

It  should  be  noted  that  the  random  phase  model  of  Equation  (1) 
represents  an  entire  ensemble  of  sea-surface  time  histories.  Here, 
ensemble  is  used  in  the  sense  of  statistical  mechanics  and  indicates  an 
assignment  of  probability  to  the  functions  of  time  which  describe  possible 
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sea-surface  histories.  Further,  any  specific  time  history  resulting  from 
this  representation  will  be  a deterministic  function,  arbitrarily  selected 
from  the  infinite  number  of  such  time  histories  that  make  up  the  ensemble. 
Detailed  theoretical  discussion  of  the  random  phase  wave  surface  model  is 
given  in  Appendices  A through  C.  There,  Equation  (1)  is  derived  from  first 
principles  and  its  relationship  to  the  more  usual  Riemann-Stielt jes 
integration  notation  is  defined. 

Equation  (1)  is  essentially  the  same  as  the  form  given  by  St.  Denis 

and  Pierson^-  as  a theoretical  ocean  wave  surface  model,  except  that  2 S(to) 

2 

is  used  instead  of  n (to)  and  the  sign  of  the  phase  angles  is  negative 
rather  than  positive.  The  St.  Denis-Pierson  paper  also  described  roughly 
how  the  integral  in  Equation  (1)  is  obtained  as  the  limit  in  quadratic 
mean  (i.i.m.)  of  a sequence  of  random  partial  sums  (see  Appendix  A for  a 
detailed  treatment).  That  is 

A 

ri ( t ) = £.i.m.  I cos  (tot— C (to) ) /2S (to) dto  (2) 

A 0 


where 


(tot— £ (to) ) 


/2S(io)dio  = £.i.m. 


Z 

k=l 


* * 
cos  ) 


A 


,2S(tok)(tok-tok_1) 


max  | iok-(0k  1 1 -+  0 

N ■*  °°  (3) 

In  the  right  hand  side  of  Equation  (3)  the  summation  is  over  all  subinter- 
vals appearing  in  the  partition 

0 = (oQ  < to1  . . . < Vl  < - A 


* 


of  the  interval  [0,A],  and  iok 


is  an  arbitrary  point  in  (tok  ,tok)  . 
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It  follows  from  the  uniform  distribution  of  phase  angles  that  the 
random  phase  model  has  a Gaussian  character.  As  shown  in  Appendix  C,  the 
model  has  the  following  statistical  properties: 


Property  1 
Property  2 


For  each  time  t,  the  expected  value  (mean  value)  of 
the  wave  surface  elevation  n(t),  is  zero. 


For  each  time  t,  the  variance  of  r|(t)  is 
The  RMS  value  (standard  deviation)  is 

1/2 


Jo 


S(w)duj. 


a = 


f 

Jo 


S(w)du) 


Property  3 - The  wave  process  is  stationary  and  has  the  covariance 
function 


R(t) 


-f 

•JO 


S(u)  cos  (u)T)du 


Property  4 - For  each  t,  n(t)  has  a normal  (Gaussian)  distribution 

2 

with  mean  zero  and  variance  a . 

Property  5 - For  distinct  times,  t^,  ....  t^,  the  joint  distribution 

of  (n(t, ) , ...,  n(t  ))  is  multivariate  normal  with 
i n 

mean  vector  0 and  covariance  matrix  M,  with  elements 

Mij  ‘ "“rV- 

Property  6 - The  process  n(t)  is  ergodic. 


These  properties  imply  that  the  random  phase  model  is  a zero  mean, 
stationary  Gaussian  process.  The  process  is  also  ergodic,  which  means 
that  its  statistical  moments  can  be  obtained  as  averages  ovei  o single 
(infinite)  time  history.  A possible  point  of  confusion  between  the  energy 
spectrum,  S (to) , and  the  more  usual  statistical  term,  spectral  density 
function,  may  be  resolved  as  follows.  Consider  a stationary  Gaussian 
process,  n(t),  with  covariance  function  R(t).  The  Fourier  transform 
fn(u)  of  R(x)  is  called  the  spectral  density  function  of  ri(t).  Thus, 
f^(w)  is  defined  over  the  interval  (-00,00) . The  wave  energy  spectrum  S(u)) 
is  then  defined  by 


S(w) 


2fn0o) 
' 0 


if  O)£[0,°°) 

elsewhere. 


(A) 
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The  preceding  discussion  summarizes  the  theoretical  features  of  the 
continuous  time,  continuous  frequency,  random  phase  model  of  ocean  waves. 
What  follows  is  a description  of  the  authors'  procedure  for  digitally 
synthesizing  ocean  wave  time  histories  from  the  theoretical  model.  It 
turns  out  that  the  abstract  definition  of  the  model  given  above  can  be 
directly  exploited  to  obtain  the  digital  results. 

In  order  to  implement  Equation  (1)  on  a digital  computer,  the  infinite 
upper  limit  of  integration  must  be  replaced  by  a finite  value  and  the 
limiting  sequence  of  partial  sums  in  Equation  (3)  must  be  replaced  by  a 
single  sum  involving  N frequency  components.  Neither  of  these  restrictions 
poses  a real  problem  in  practice.  Ocean  waves  are  basically  low  frequency 
phenomena  and  contain  very  little  energy  at  frequencies  beyond  tt  rad/sec. 
Also,  the  wave  spectrum  can  be  approximated  as  closely  as  desired  by  the 
use  of  a finite  number  of  frequencies.  Thus,  by  choosing  a finite  maximum 
frequency  Wc,  and  a large  but  finite  number  of  frequency  subintervals,  we 
obtain  the  following  approximation  (to  Equation  (1)) 

N 

n(t)  = Xcos  (w£t-e(u>*))  >/2S(w*)  (wk-wk_i)  (5) 

k=l 


with 


0 = < a),  . . . < U)  , < Cl)  = co 

0 1 N-l  N c 

* 

By  definition,  eacli  can  be  chosen  arbitrarily  from  the  subinterval 

* 

and  each  must  be  a random  phase  angle  having  an  uniform 

distribution  on  [0,2tt).  Standard  algorithms  exist  for  generating  random 
numbers  on  digital  computers.  Discussion  of  such  algorithms  with  partic- 
ular emphasis  on  the  quality  of  random  numbers  generated  by  the  CDC  6700 
computer  system  is  given  in  Appendix  D.  However,  by  using  some  such 
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generator,  the  can  be  chosen  at  random  from  the  interval 
the  e(iok)  and  be  chosen  from  an  uniform  distribution  on  [0,2ir). 

For  convenience,  the  in  Equation  (5)  can  be  equally  spaced, 
the  approximate  wave  system  becomes 


and 

Then 


1/2 


n< 


v—*  * ■*■/  * * * 

t)  = y [2S(uk)Acj]  cos  (wkt-e(o)k)) 


k=l 


(6) 


where 


Aco  = wk  - o>k  ^ , k = 1 , . . . N 


(7) 


If  we  denote 


* 1'2 
Ck  = [2S(wk)Au>] 


(8) 


and 


e(<V  " ek 

Equation  (6)  can  be  rewritten  in  the  form 

N 

n(t)  = 2\  cos  0\t_ek)  (9) 

k=l 

Equation  (9)  is  then  used  to  digitally  generate  time  histories  in  the 
computer  program.  Although  Equation  (9)  is  not  a Fourier  series,  it  does 
contain  only  a finite  number  of  periodic  components  and  thus,  in  prinicple, 
it  will  be  periodic.  However,  because  the  frequency  components  are 


1 


A 
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selected  randomly  within  each  subinterval,  the  probability  is  zero  that 
two  distinct  components  will  be  harmonics.  Thus,  the  wave  period  resulting 
from  the  use  of  such  components  will  be  determined  by  the  number  of 
components  employed  and  the  number  of  significant  digits  used  to  represent 
each  component's  period.  Based  on  the  fifteen  significant  digits  used  in 
CDC  6700  computations,  a simulation  based  on  N wave  frequency  components 
over  the  interval  [0 ,tt) , each  having  a period  containing  significant 
digits,  would  have  a period 


T = 


N 

X 

k=l 


10 


15N 


sec 


(10) 


Thus,  assuming  the  selection  of  a moderate  number  of  frequency  components, 
for  example  50,  the  wave  period  would  be  roughly 


T = 10750  sec 


= 3.17  x 10742  yr 


(ID 


Hence,  for  all  practical  purposes,  the  wave  period  is  infinite  so  that  the 
wave  system  may  be  considered  to  be  nonperiodic.  The  advantage  of  the  wave 
energy  integral  representation  becomes  clear  upon  noting  that  a Fourier 
series  model  containing  100  frequency  components  over  [ 0 , tt)  would  have  a 
pe  iod  of  only  200  sec. 


DESCRIPTION  OF  SIMULATED  OCEAN  WAVE  RECORDS 
The  idealized  character  of  a continuous  parameter  stochastic  process 
cannot  be  realized  in  a digital  simulation.  The  simulated  wave  time 
histories  must  be  based  on  a discrete  approximation  to  the  theoretical 
wave  spectrum.  In  addition,  only  a finite  number  of  wave  frequency 
components  and  random  phases  can  be  employed.  It  is  thus  essential  to 
investigate  how  well  the  simulated  records  approximate  the  theoretical 
ones.  In  this  respect,  the  effects  of  both  approximation  error  and 
sampling  variability  are  significant. 


r 1 

■ 

I 1 
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The  computer  program  developed  herein  allows  the  user  to  input  the 
wave  spectrum  model  of  his  choice.  It  is  merely  required  that  the 
spectral  ordinate  be  input  at  each  discrete  frequency  of  interest.  Hence, 
the  model  can  be  easily  adopted  to  model  wave  systems  exhibiting 
empirical  spectra  such  as  the  station  India  spectra  used  in  many  ship  per- 
formance studies.  Likewise,  the  model  readily  admits  theoretical  spectra 
such  as  the  International  Towing  Tank  Conference  (ITTC)  recommended  two- 

parameter  Bretschneider  spectrum.  In  the  numerical  examples  that  follow, 

7 

a 30-knot  wind  spectrum  developed  by  Pierson  and  Moskowitz  was  used  to 
generate  the  required  wave  processes.  The  frequency  range  selected  was 
[0.2, 2. 2]  rad/sec.  The  Pierson-Moskowitz  spectrum  is  defined  in  units  of 
ft^-sec  by 


S(w) 


0.00810  e-0.749/(VK-oj)4 
u>5 


(12) 


where  u)  = frequency  (rad/sec) 

VK  = wind  speed  (knots) 

The  spectrum  is  a single  amplitude  spectrum  and  is  continuous  over  positive 
frequencies.  As  implemented  in  this  report,  the  area  under  the  spectrum 
will  equal  the  variance  of  the  generated  wave  process.  Figure  1 shows  the 
energy  spectrum  based  on  the  Pierson-Moskowitz  formula  for  a wind  speed  of 
30  knots.  It  should  be  noted  that  any  desired  spectrum  can  be  used  in  the 
simulation.  The  Pierson-Moskowitz  spectrum  was  selected  only  for 
convenience. 

In  order  to  compute  wave  surface  time  histories  from  Equation  (9) , the 
number  of  frequency  components  N must  be  selected  in  a manner  that  is 
dependent  upon  the  intended  application.  In  many  instances  where  strict 
adherence  to  the  normal  distribution  of  amplitudes  is  important,  the 
minimum  number  of  frequency  components  should  be  sufficiently  large  for 
the  central  limit  theorem  to  hold.  The  normality  property  was  investigated 
by  applying  chi-square  goodness-of-f it  tests  to  samples  of  wave  amplitude 
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Figure  1 - Pierson-Moskowitz  Wave  Spectrum  at  30  Knots 


points  corresponding  to  a number  of  frequency  components  N.  Since  inde- 
pendent sampled  points  are  required  for  application  of  the  usual  chi- 
square  tests,  such  tests  cannot  be  performed  straightforwardly  on 
correlated  wave  histories.  Independent  wave  amplitudes  were  obtained  by 
sampling  the  wave  ensemble  600  times  at  the  fixed  time  t = 0.  This  pro- 
cedure was  applied  to  each  of  ten  samples  consisting  of  5,  10,  15,  25,  50, 
and  100  wave  frequency  components,  respectively.  Thus,  ten  chi-square 
tests  were  performed  for  each  set  of  frequency  components.  The  number  of 

class  intervals  was  selected  to  maximize  the  power  of  the  tests  in 

8 9 

accordance  with  the  procedure  of  Mann  and  Wald  as  detailed  by  Williams. 
Following  this  procedure,  the  number  K,  of  class  intervals  depends  upon  the 
number  of  points  N and  the  significance  level  a,  according  to  the  formula 


K = b 


(SV1) 


where  N is  the  number  of  sample  points  and  a is  defined  by 


no 

it175  / 


_ 1 
2 , 

e dy  = a 


Fifty-two  class  intervals  were  used  and  the  data  were  grouped  into  52 

equiprobability  classes  which  yielded  51  degrees  of  freedom.  The  results 

from  the  chi-square  tests  are  shown  in  Table  1.  Based  on  the  significance 

level  a = 0.10,  no  more  than  one  sample  out  of  ten  would  be  expected  to 

2 

exceed  the  critical  value,  x = 64.295,  if  the  data  were  representative  ot 

2 

a normal  distribution  with  mean  zero  and  variance  a . For  five  frequency 
components,  all  ten  of  the  computed  values  exceed  the  critical  value.  For 
twenty-five  or  more  components,  none  of  the  ten  chi-square  values  exceeds 
the  critical  value.  Chi-square  tests  were  also  computed  on  the  data  by 
using  the  fifty-two  class  intervals  but  without  grouping  the  data  into 
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equally  probable  intervals  and  the  same  trends  were  observed.  Hence,  as 
few  as  twenty-five  frequency  components  are  sufficient  to  guarantee  con- 
vergence of  the  random  phase  model  to  the  normal  law.  However,  this 
conclusion  is  strictly  valid  only  for  the  range  of  variation  observed  in 
the  samples,  namely  ± 3o.  More  frequency  components  may  be  required  to 
imply  normality  for  the  tails  of  the  sample  distributions. 

An  example  of  a generated  ocean  wave  time  history  is  shown  in  Figure  2. 
This  time  history  is  a 30  min  record  based  on  100  frequency  components 
over  the  frequency  range  [0.2, 2. 2]  rad/sec.  A sampling  rate  of  1 sample 
per  second  was  used.  Figure  3 shows  a sample  probability  density  function 
of  the  wave  amplitudes  from  the  30  min  record.  Close  agreement  is  ob- 
served between  the  computed  wave  data  and  the  theoretical  normal  distribu- 

2 

tion  based  on  mean  zero  and  variance  0 . The  usual  chi-square  goodness-of- 
fit  test  could  not  be  performed  on  the  time  history  data  because  adjacent 
sample  points  are  not  independent. 

An  important  property  of  the  wave  record  is  its  mean  square  value  or 
variance.  The  theoretical  variance  of  the  generated  wave  process  is  given 
by  the  area  under  the  wave  spectrum.  That  is 


o 


2 


S(w)dw 


(15) 


where  a and  b are  the  minimum  and  maximum  frequencies,  respectively,  in 
rad/sec,  and  S(w)  is  the  energy  spectrum.  The  variance  of  the  theoretical 
wave  process  based  on  the  continuous  form  of  the  Pierson-Moskowitz 
spectrum,  where  b * 2.2  and  a = 0.2,  is  17.287.  The  variance  of  the 
generated  wave  process  was  computed  as  a function  of  record  length  and 
number  of  frequency  components.  Figure  4 shows  the  computed  variance  of 
additional  wave  records  based  on  25,  50,  and  100  frequency  components, 
respectively,  and  based  on  time  histories  between  0 and  30  min  in  length. 
The  horizontal  line  represents  the  theoretical  variance  based  on  the  30- 
knot  Pierson-Moskowitz  spectrum.  In  time  histories  of  length  less  than 
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Figure  4 - Variance  of  Simulated  Wave  Records 


35 


16 


five  min,  the  computed  variance  fluctuates  considerably  and  does  not 
approximate  the  theoretical  variance  too  closely.  As  the  length  of  record 
increases,  the  approximated  quality  improves  and  the  error  is  within 
± 5 percent  for  records  no  shorter  than  10  min.  As  expected,  the  variance 
approximation  improves  as  the  number  of  frequency  components  is  increased. 
However,  25  frequency  components  are  seen  to  yield  reasonable  accuracy. 

The  theoretical  covariance  function  of  the  generated  wave  process  is 
given  by  (Appendix  C) 


R(x ) = J (S(w)  cos  (uyt))dT  (16) 

0 

This  theoretical  covariance  function  is  shown  in  Figure  5 along  with  the 
sample  covariance  function  computed  from  the  30  min  record  (100  frequen- 
cies). The  two  curves  are  in  very  close  agreement  for  lags  up  to  60  sec. 

Figure  6 shows  a comparison  of  a sample  spectrum  computed  from  the 
30  min  record  and  the  theoretical  Pier son-Mo skowitz  spectrum.  Close 
agreement  is  observed  between  the  two  curves,  although  the  sample  spectrum 
shows  some  oscillation  due  to  sampling  variability. 

THE  COMPUTER  PROGRAM 

PROGRAM  DESCRIPTION 

This  section  describes  the  computer  program  developed  for  generating 
and  plotting  the  random  phase  wave  records.  WAVESIM,  the  main  program, 
functions  as  an  executive  program.  It  reads  input  parameters  that  must  be 
supplied  by  the  user.  Table  2 explains  the  parameters  used  in  the  pro- 
gram. The  subroutines  that  follow  are  called  from  the  main  program,  and 
the  input  data  is  transferred  to  these  subroutines  for  computation.  A 
listing  of  the  program  appears  in  Appendix  E.  The  program  consists  of  the 
subroutines  SPECTRA,  WAVES,  ARAMEAN,  VAR,  and  PLOT  W,  which  are  discussed. 
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FREQUENCY  Irad/wc) 

Figure  6 - Estimated  Spectral  Density  Function  (30  Minute  Record) 
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TABLE  2 - COMPUTER  PROGRAM  NOMENCLATURE 
Number  of  cases 

Number  of  nonzero  frequency  subdivisions 
M + 1 

Time  increment  (sec) 

Minimum  frequency  (rad/sec) 

Maximum  frequency  (rad/sec) 

Length  of  the  wave  process  (min) 

Use  1 if  plot  is  desired.  Use  2 if  not 

Frequency  increment  (rad/sec),  DELW  = (OMEGMAX-OMEGMIN) /M 
Length  of  the  wave  process  (sec) 

Total  number  of  points  computed  for  wave  process, 

NO  = TINSEC/DELT  + 1 

Frequency  subdivision 

Actual  frequency 

Pierson-Moskowitz  spectral  ordinate  based  on  actual 
frequency 

Pierson-Moskowitz  spectral  ordinate  based  on  frequency 
subdivision 

Random  phase  angle  from  uniform  distribution  on 
[0,2tt1  (rad) 

Time  (sec) 

Wave  Process  as  a function  of  time 

Time  average  of  the  wave  process 

Variance  of  the  computed  wave  process 

Area  based  on  the  continuous  form  of  the  Pierson- 
Moskowitz  spectrum 

Wind  speed  (knots) 
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SPECTRA:  This  subroutine  calculates  the  spectral  density  function 

using  the  Pierson-Moskowitz  formula  given  by: 


x 0.00810  -0. 749/  (VK'oj) 

S(co)  = ^ e 

U) 


where  S(io)  = ordinate  of  spectral  density 
u)  = circular  frequency 

VK  * wind  speed  (knots) 

This  formula  requires  a value  for  wind  speed  in  knots  (VK)  which  is  read 
as  input  in  the  program.  The  frequency  subdivisions  and  the  actual  fre- 
quencies are  generated  in  this  subroutine.  The  spectral  density  function 
is  defined  only  over  non-negative  frequencies.  This  subroutine  is 
optional  as  the  user  may  substitute  any  spectrum  S(w)  that  is  desired. 

WAVES:  This  subroutine  computes  (using  Equation  (9))  the  random  wave 

process  using  the  formula: 


n(t)  = 


V”*  * 1/2  * 

/ [2S(to^)Aco]  cos  (w^t-Cj^) 


I 


where  S(w)  = the  spectral  ordinate  values 
Aw  = frequency  increment  (rad/sec) 

•k 

0)^  = actual  frequency  (rad/sec) 

t = time  (sec) 

= random  phase  angle 

N = number  of  discrete  frequency  components  to 
be  employed 

The  values  for  time  in  seconds  and  the  random  phase  angles  are  generated 
within  this  subroutine. 


ARAMEAN:  This  subroutine  calculates  the  mean  value  of  the  generated 

wave  process.  The  area  under  the  continuous  form  of  the  Pierson- 
Moskowitz  spectrum  is  also  computed. 


VAR:  This  subroutine  computes  the  variance  of  the  simulated 

time  history. 

PLOT  W:  This  subroutine  plots  the  generated  wave  on  the  CALCOMP  936 
plotter.  The  value  of  IFPLOT  (input  to  the  main  program)  determines 
whether  or  not  this  routine  is  used.  If  a plot  is  desired,  the  program 
generates  a magnetic  tape  from  which  plots  can  be  obtained  from  the 
CALCOMP  plotter. 

PREPARATION  OF  INPUT 

The  input  for  calculating  and  plotting  a wave  history  consists  of 
three  cards  which  are  described  in  Table  3.  There  are  seven  parameters, 

N,  M,  MM,  DELT,  OMEGMIN,  OMEGMAX,  and  TINMIN,  on  the  first  card.  The 
parameter  on  the  second  card  is  IFPLOT.  All  integer  names  must  be  right 
justified  since  blank  columns  are  interpreted  as  zeros.  Floating  point 
names  must  include  a decimal  point.  If  the  Pierson-Moskowitz  spectrum  is 
used,  the  third  card  nust  contain  the  wind  speed  parameter  VK.  Table  4 
shows  an  example  of  the  input  with  and  without  the  plot  routine.  If 
another  case  follows,  an  additional  set  of  data  (Cards  1-3)  is  needed. 

The  last  case  should  have  a value  of  1 for  N,  in  order  to  terminate  the 
program. 

EXPLANATION  OF  OUTPUT 

The  output  consists  of  eight  sections  of  computed  information.  A 
sample  printout  is  shown  in  the  following  sections.  Section  1 consists  of 
the  input  parameters.  Section  2 shows  the  frequency  subdivisions,  actual 
frequencies,  and  spectral  ordinate  values  of  the  Pierson-Moskowitz 
spectrum  based  on  the  actual  frequencies  between  0.2  and  1.0.  The  random 
phase  angles  are  listed  in  Section  3.  Section  4 shows  the  first  136 
points  of  the  generated  wave  process  as  a function  of  time  in  seconds.  In 


22 


I 


TABLE  3 - INPUT  FOR  MAIN  PROGRAM 


Card 

Columns 

Format 

FORTRAN 

Designation 

Definition 

1 

1-  5 

15 

N 

Number  of  cases 

6-10 

15 

M 

Number  of  frequency 
subdivisions 

11-15 

15 

MM 

MM  = M + 1 

16-25 

F10. 2 

DELT 

Time  increment  (sec) 

26-35 

F10. 2 

OMEGMIN 

Minimum  frequency  (rad/sec) 

36-45 

F10.2 

OMEGMAX 

Maximum  frequency  (rad/sec) 

46-55 

F10.2 

TINMIN 

Length  of  the  wave 
process  (min) 

2 

1-  5 

15 

IFPLOT 

Use  1 if  plot  routine  is 
desired.  Use  2 if  not. 

3 

1-10 

F10. 2 

VK 

Wind  speed  (knots)  for  the 
Pierson-Moskowitz  formula 
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SECTION  1 


NUMBER  OF  INPUT  CASES  1 

TOTAL  NUMBER  OF  POINTS  COMPUTED  FOR  WAVE  PRUCFSS  1801 
NUMBER  OF  FREQUENCIES  USED  101 

MINIMUM  FREQUENCY  IN  RAO  I ANS/SECOND  .300 

MAXIMUM  FREQUENCY  IN  RADIANS/SECOND  2.300 

FREQUENCY  INCREMENT  IN  RADIANS/SECOND  .020 

TIME.  INCREMENT  IN  SECONDS  1.000 

LENGTH  OF  wave  PROCFSS  IN  MlNUTtS  30.000 

LENGTH  OE  WAVE  PROCESS  IN  SECONDS  1800.000 

WIND  SPEED  IN  KNOTS  30.000 


2S 


SECTION  2 


FREQUENCY  ACTUAL  SPECTRAL  ORDINATE  BASED 

SUBDIVISION  FREQUENCY  ON  ACTUAL  FREQUENCY 


poo 

.200 

220 

.221 

240 

.252 

260 

.261 

P80 

.294 

300 

.302 

320 

.337 

340 

.355 

360 

.376 

380 

.382 

400 

.405 

420 

.429 

440 

.449 

460 

.461 

480 

.495 

500 

.515 

520 

.538 

540 

.555 

.000 

.000 

.000 

.000 

.000 

.00? 

.161 

.737 

2. 729 

3.619 

8.729 
16.415 
23.645 
27.778 
37.897 
41.606 
44.061 
44.648 


.560  .569  44.478 
.580  .595  42.975 
.600  .602  42.340 


.620 

.624 

• 640 

.652 

.660 

.664 

.680 

.688 

.700 

.711 

.720 

.735 

.740 

.743 

.760 

.779 

.780 

.784 

.800 

.807 

.820 

.836 

.840 

.858 

.860 

.862 

.880 

.893 

.900 

.917 

.920 

.931 

.940 

.952 

.960 

.974 

.980 

.983 

1.000 

1.001 

40.052 

36.461 
34.920 
31.710 
28.772 
25.884 
24.945 
21.077 
20.541 
18.463 
16.048 

14.461 
14.186 
12.221 
10.883 
10.222 

9.262 

8.384 

8.021 

7.396 
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SECTION  3 


RANDOM  PHASE  ANGLES 


4. 296128 

5.648159 

2.762822 

.761861 

5.094957 

4. 373786 

i. 421783 

2. 834857 

4.681262 

5.759672 

2. 743953 

4.526702 

4.427508 

5.621661 

2.713662 

4. 247852 

4.246169 

4.790832 

2.210933 

1.973443 

2. 920811 

.006001 

5.817229 

5.293417 

2.636457 

4.  952821 

5.108381 

4.354680 

4.912478 

4.862324 

5. 426944 

.579937 

2.156256 

2.009325 

1.669190 

2. 060736 

2.832687 

2.521395 

3.426094 

.471848 

5. 234771 

4.813573 

3.441344 

1.295194 

3.138264 

5. 170657 

6.179425 

6.253782 

2.337308 

6.022332 

3. 169101 

3.479571 

3.249147 

.953576 

3.601916 

4.  185023 

6.236676 

.838576 

1.062601 

5.199075 

3. 002502 

5.552497 

2.112288 

5.  060934 

. 089845 

3. 887805 

2.693324 

4 • 536166 

1.824741 

5.764042 

. 017589 

2.827710 

5.295909 

1.889345 

1.405937 

4. 277415 

.783155 

.708669 

5.368163 

2. 675637 

2. 326224 

2.174309 

3.432989 

3. 190898 

1.607190 

1. 256245 

4.906084 

1.297543 

.872652 

2.645112 

2. 370641 

1.055278 

4.611419 

1.917550 

2.652305 

3. 891230 

3. 942509 

5.905635 

4.989108 

i. 775157 

4.563567 

SECTION  4 
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Section  5 the  mean  and  variance  of  the  wave  process,  the  area  based  on  the 
continuous  form  of  the  Pierson-Moskowitz  spectrum,  and  the  minimum  and 
maximum  values  of  the  wave  process  are  printed. 

OPERATING  COMPUTER  INFORMATION 

The  WAVES IM  program  is  coded  in  FORTRAN  IV  and  has  been  used  on  the 
CDC  6700  computer  system.  However,  the  program  should  be  easily  adapted 
to  any  other  system  provided  an  uniform  random  number  generator  is  avail- 
able. Table  5 shows  the  sample  deck  structure  for  the  CDC  6700  system 
with  and  without  the  plot  routine.  The  number  of  words  of  computer  memory 
and  the  computer  time  required  for  the  program  will  depend  on  the  number 
of  frequencies  used  and  the  length  of  the  wave  process.  For  example,  a 
typical  run  involving  the  computation  and  plotting  of  a 30  min  record 
based  on  100  frequency  components  and  sampled  at  1 sample/sec,  required 
66.5  system  sec  on  the  CDC  6700. 
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TABLE  5 - SAMPLE  DECK  STRUCTURE  FOR  THE  CDC  6700  SYSTEM 


1.  Compilation  and  Execution  with  PLOT  ROUTINE,  IFPLOT  “ 1. 
JOBNAME,CM60000,MT1,P3,T100.  *USERS  JOB  CARD* 

CHARGE, CXXX.XXXXXXXXXX.  *USERS  CHARGE  CARD* , *XXXXXXXXXX  IS 

VSN,  TAPEN,  TAPENAME=*SLOT  . USERS  JOB  NUMBER*,  *CXXX  IS  USERS  ID* 

REQUEST , TAPEN , HI . (TAPENAME , SLOT /RING) 

FTN(T) 

ATTACH, CALC 9 36. 

LDSET (LIB»CALC936) 

LGO. 

RETURN (TAPEN) 

7/8/9  END  OF  RECORD 
FORTRAN  DECK 
7/8/9  END  OF  RECORD 
DATA 

6/7/8/9  END  OF  FILE 


2.  Compilation  and  Execution  without  PLOT  ROUTINE,  IFPLOT  “ 2. 


JOBNAME,CM60000,P3,T100. 
CHARGE , CXXX , XXXXXXXXXX 
FTN(T) 

LGO. 

7/8/9  END  OF  RECORD 
FORTRAN  DECK 
7/8/9  END  OF  RECORD 
DATA 

6/ 7/8/9  END  OF  FILE 
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APPENDIX  A 

RANDOM  FUNCTION  THEORY 


A few  basic  facts  concerning  random  functions  are  required  to  develop 
and  interpret  the  random  phase  ocean  wave  model.  A brief  sumi^ary  of  this 
theory  is  given  below  for  immediate  reference.  More  rigorous  discussion 
and  proofs  may  be  found  in  standard  textbooks  such  as  Doob^  or  Lukacs.^ 

A random  function  X is  a function  whose  domain  is  a parameter  set  T, 
and  whose  range  is  a set  of  random  variables,  possibly  complex  valued. 
Thus,  a random  function  may  be  considered  as  an  indexed  set  of  random 
variables,  {Xt|t£T}. 

In  practice,  the  set  T will  usually  consist  of  discrete  or  continuous 
subsets  of  the  real  line.  When  only  continuous  subsets  are  considered,  it 
is  conventional  to  write  each  member  as  X(t)  and  to  denote  the  process  by 
{X(t) |teT}. 

For  each  t,  the  random  variable  Xt  is  determined  by  its  probability 
distribution  function,  Fxt(x),  as  in  elementary  probability  theory.  For 
each  Xt,  moments  of  arbitrary  order  can  be  computed,  whenever  they  exist, 
by 


E X"  = J Xn 


dFxt(x) 


(A.  1) 


where  E denotes  the  mathematical  expectation  or  mean  value  operator.  At 
any  two  different  times  t^  or  the  covariance  function  between  the  twc 
random  variables  is  defined,  whenever  it  exists,  by 


COV 


(V\)  ‘"[(VS)  (V\)] 


(A. 2) 


A R(t1,t2) 
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An  important  class  of  random  functions,  useful  in  both  theoretical 

work  and  applications,  is  characterized  by  the  fact  that  the  principal 

properties  of  the  indexed  random  variables  are  invariant  with  respect  to 

the  location  of  the  random  variables  within  the  indexed  set.  Such 

processes  are  called  stationary.  A process  {Xt|teT}  is  said  to  be  wide 

sense  stationary  if 

2 2 

1.  VteT  EX^,  EXt  exist  and  are  independent  of  t. 

2.  Vt^,t2£T,  the  covariance  function  R(t^,t2)  exist  and  are 

invariant  under  any  translation  of  the  t axis. 

2 

Processes  for  which  EXt  exists  are  called  second-order  processes.  It 

follows,  by  Schwartz's  inequality,  that  for  any  second-order  process, 

X , COV  /X  , X \ < °°.  It  is  clear  from  the  second  property  of  stationary 
V C1  C2  / 

processes  that  the  covariance  function  can  depend  only  on  the  time 
difference  between  t^  and  t2  and  not  on  t^  or  t2  separately.  For  this 
reason,  it  is  customary  to  denote  T = t2  - t^,  and  write 

R(T)  - COV 

= COV  ^Xt  ,Xt  +T\  (A. 3) 

The  notion  of  integration  for  random  functions  can  be  defined.  How- 
ever, as  in  the  usual  Riemann-Stieltjes  integral,  the  concepts  of  con- 
vergence to  a limit  point  and  continuity  are  needed.  Let  {X^}  be  a 
sequence  of  second  order  random  variables.  A random  variable  X is  said  to 
be  the  limit  in  quadratic  mean  of  {X^}  if 

1.  X is  second  order 

2.  Z. i.m.  E| (X  -X) |2  = 0 
n -*■  00 

We  write 

X = l. i.m.  X (A. 4) 

n 00  n 

and  we  say  that  X^  converges  in  quadratic  mean  to  X. 
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Let  {X(t)|teT}  be  a second  order  continuous  parameter  stationary 
process.  Then  the  process  is  said  to  be  continuous  on  T if 

VteT  Jl.i.m.  X(t+x)  = X(t)  (A. 5) 

T ■>  0 

A second  order  stationary  process  is  continuous  in  the  sense  of  Equation 
(A. 5)  if  and  only  if  its  covariance  function  is  continuous  in  the  usual 
sense  at  T = 0.  Continuous  processes  have  the  property  that  they  do  not 
fluctuate  too  wildly  and  are  thus  amenable  to  integration.  The  stochastic 
integral  may  be  defined  in  a fashion  analogous  to  the  Riemann-Stieltjes 
integral,  the  major  difference  being  that  the  real  line  convergence  of  a 
sequence  of  real  sums  must  be  replaced  by  some  mode  of  stochastic  con- 
vergence for  the  corresponding  sequence  of  random  sums.  The  following 
definition  can  be  given  based  on  quadratic  mean  convergence. 

Let  X(t)  and  Z(t)  be  any  two  second  order  continuous  processes  defined 
on  the  finite  interval  [A,B] . Consider  the  sequence  of  partitions  of 
[A,B]  given  by 


A 

n 


n,  1 


< t 

n,n 


B 


Then,  for  every  n,  let 


S = 
n 


jrX<k)1Z(tA,k>  - 


k-1 


where 


Assume  further  that 


t . e[t  , , ,t  . ] 
n,k  n,k-l’  n,k 
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ft.i.m.  max  t . - t . , = 0 

n co  k n’k  n»k"1l 


Then,  if  there  exists  an  unique  random  variable  S,  such  that 

S = £.i.m.  S 

n 

n 00 

we  call  S the  integral  over  [A,B]  of  X(t)  with  respect  to  Z(t)  and  denote 
it  by 


.B 

S = | X(t)  dZ( t) 

'A 


f 


(A.  6) 


It  can  be  shown  that  if  X(t)  and  Z(t)  are  independent  VteT,  the  integral 
exists  if  and  only  if 


B B 


ESI 


= J J Rz(s,t)Rx(s,t)dsdt 


(A.  7) 


A A 


where  Rz  and  Rx  are  the  respective  covariance  functions. 

Special  cases  of  the  integral  occur  when  either  X(t)  or  Z(t) 
degenerates  into  a deterministic  function.  The  integral  can  be  defined 
for  an  infinite  interval  by  the  standard  limiting  procedure  used  for 
Riemann-S t iel t j es  integrals.  The  stochastic  integral  is  linear  and  homo- 
geneous in  the  integrand  and  additive  in  the  domain  of  integration. 

It  is  often  desirable  to  perform  integrations  with  respect  to 
processes  possessing  orthogonal  increments.  A process  {Xt|teT}  is  said  to 
have  orthogonal  increments  if 


(a)  E | Xt-Xg  | 2 < °°  Vs.teT 


(A. 8) 
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do 


(b)  If  the  two  intervals  [t^.t^,  [t^t^] 

not  overlap,  then  the  increments  /Xt  -X  \ 
^Xfc  -Xt  j are  orthogonal,  i.e.,  ^ ^ 


and 


(VS)  (VV 


(A.  8) 

(cont . ) 


where  the  bar  denotes  complex  conjugation.  The  importance  of  processes 
possessing  orthogonal  increments  is  due  largely  to  the  celebrated  Spectral 
Representation  Theorem.  This  theorem  states  that  any  wide  sense  station- 
ary, continuous  stochastic  process  Y(t)  can  be  represented  as  the 
Fourier-Stieltjes  integral 


Y(t) 


OO 

J eiUt  dX(co) 

— OO 


(A.  9) 


where  X(t)  is  a process  with  orthogonal  increments  and  E | dX(aj)  | = dG(to)  . 

Further,  X(t)  when  properly  normalized,  is  uniquely  specified  by  this 
relationship. 

The  order  of  expectation  and  integration  can  be  interchanged  in  the 
stochastic  integral  so  that 


EY(t)  = 


d (EX(w) ) 


(A. 10) 


whenever  the  right  hand  side  of  Equation  (A. 10)  exists.  Since  X(oj)  is  a 
random  variable  for  each  co,  Y(t)  consists  of  an  infinite  number  of  time 
histories  over  the  range  (-00,00) , and  each  time  history  is  defined  by  a set 
of  simultaneous  outcomes  of  the  random  variables  (X(oj)  |a)e(-°0,00) } . Such 
time  histories  are  called  sample  functions  or  realizations  of  the  process 
Y(t).  If  Y(t)  is  a continuous  (ft.i.m.)  process,  each  sample  function  will 
be  a deterministic  function  and  will  be  continuous  except  possibly  on  sets 
of  zero  probability. 
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A stochastic  derivative  can  also  be  defined  for  a second-order 
process,  Y(t).  Analogous  to  the  usual  derivative  definition,  we  define 


Y'(t)  = Z. i.m.  Y(t+T)  ~ Y(t) 

T 

T ->  0 


whenever  the  limit  exists.  If  R(s,t)  is  the  covariance  function  of  Y, 
Y(t)  exists  whenever 


92R(s,t) 

3s9t 


■] 


s=t 


exists.  The  function  Y'(t)  is  called  the  quadratic  mean  derivative  of 
Y (t) . 
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APPENDIX  B 

STOCHASTIC  INTEGRAL  REPRESENTATION  OF  THE  RANDOM 
PHASE  MODEL 


The  mathematical  form  of  the  random  phase  wave  surface  model  proposed 
by  St.  Denis  and  Pierson  may  be  unfamiliar  to  the  reader.  The  model  in- 
volves the  representation  of  a random  function  by  what  appears  to  be  an 
usual  integration  process,  but  which  contains  a square  root  symbol  over 
the  differential  frequency  element.  Appendix  A contains  the  background 
material  required  to  properly  interpret  the  integral  of  a stochastic 
function.  In  this  appendix,  the  spectral  representation  theorem  is  em- 
ployed to  derive  the  random  phase  model  in  terms  of  a random  Fourier- 
Stieltjes  integral  which,  notationally , does  not  require  use  of  the  square 
root  symbol.  However,  the  square  root  remains  implicitly  and  its 
importance  to  the  representation  is  shown.  It  will  be  shown  that 
physically,  the  square  root  implies  that  the  wave  amplitude  is  proportional 
to  the  square  root  of  cumulative  wave  energy. 

Assume  that  the  sea  surface  is  in  a stationary  state  and  we  desire  to 
represent  the  wave  surface  elevation  r] ( t ) as  a stationary  Gaussian  process. 
According  to  the  spectral  representation  theorem,  it  must  be  possible  to 
represent  the  real  process  n(t)  in  the  form 


n(t) 


Re 


J eia)t  dZ(to) 


(B.l) 


where  Z(to)  is  a process  having  orthogonal  increments.  We  shall  construct 
Z(to)  as  follows.  Let  f (to) , -°°  < to  < °°,  denote  the  two-sided  wave  energy 

i _iekl 

spectral  density  function.  Consider  the  sequence  (e  ( where  each  is 
a random  variable  uniformly  distributed  over  [0,2 tt) , i.e., 


P 

r 


[o±V27T]  = h 


(B.  2) 


1 
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Define  the  random  set  function  Z on  intervals  of  the  real  line  by 

Zfa.e)  « E.i.m.^y^e-1'*  [fn(u£)(wk-U)k_i)J  (B-3) 

u^e  (a,  6) 

maxlCJk_wk_1l  0 

where  the  w.  comprise  a sequence  of  partitions  of  the  frequency  axis  and 
* * 

ioke(a,6].  Then  define  the  random  point  function  Z(u))  by 

Z(oj)  = Z^-oo.u)])  (B.4) 


and  denote 


Z(a))  = 


I 


e'ie(a)  ^n('60da 


(B.5) 


Note  that  Z(uj)  exists  for  every  to  because 


e|Z(w)  I = EZ(u))  Z(u>) 


0) 

■I 


f^(a)dot  < °° 


It  follows  easily  that  Z(oj)  has  orthogonal  increments  and  that 
E | dZ(o))  | = dG(u),  where  G(u))  is  the  spectral  distribution  function, 

defined  by 


G(w) 


L 


f^CoOda 


(B.6) 


AO 


Now,  although  Z(w)  is  not  differentiable,  because  of  the  manner  in 
which  Z(uj)  was  constructed,  we  can  write 

dZ(tu)  = e_i£(a))  /f  (u)du  (B.7) 

Thus,  the  complex  form  of  the  random  phase  model  is 


n(t)  * 


dZ(co) 


(B.8) 


with  Z(co)  given  by  Equation  (B.5).  The  real  random  phase  model  is  then 
defined  by 


n(t)  = Re 


oo 


— OO 


dZ(u)) 


(B.9) 


Equations  (B.8)  and  (B.9)  formally  represent  the  random  phase  model  in  the 
usual  manner  for  Fourier-Stielt j es  integrals.  However,  the  square  root 
remains  an  essential  part  of  Z(uj)  as  is  seen  from  Equations  (B.5)  and 
(B.7).  Equation  (B.9)  then  yields  the  integral  form  used  in  the  simulation 
program.  Because 


n(t) 


dZ(co) 


by  (B.9) 


OO 

= Re  J eia)t  e~ie(u))  by  (B.6) 

—OO 
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(o>t-e  ((jo)  ) 


/f^(ti))dco 


by  definition 


n(t)  = 2 


00 

J cos  (wt-e(a)))  /f^(ui)daj 
0 


by  symmetry. 


(idt-e(u)))  /2S(o))d(jj 


by  definition  of  S(w)  as  the  one-sided  wave  spectral  density. 

The  necessity  of  the  square  root  also  follows  from  the  spectral 
representation  theorem  which  requires 


dG^(uj)  = cjdZ(w)  )2 


(B. 10) 


where  G^u))  is  the  wave  process  spectral  distribution  function,  called 
"cumulative  energy  density"  in  the  St.  Denis-Pierson  paper.  in  most 
practical  applications,  G^(u)  is  absolutely  continuous  so  that 


dG  (ui)  = f ((jj)du) 

n n 


(B.ll) 


Thus , 


E | dZ(to)  | 2 = f (w)du) 


(B. 12) 


which  implies  that  dZ(u)),  and  hence  n(t),  in  Equation  (B.l)  must  be  pro- 

1/2 

portional  to  (dca)  . Physically,  this  result  means  that  the  squared  wave 
amplitude  must  be  proportional  to  the  cumulative  wave  energy. 
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APPENDIX  C 

STATISTICAL  PROPERTIES  OF  THE  RANDOM  PHASE  MODEL 


We  wish  to  show  that  the  random  phase  wave  surface  elevation  model 
represents  a stationary  Gaussian  process  having  zero  mean.  Let  n(t) 
denote  the  time  varying  wave  surface  elevation  and  let  S(oj)  denote  the  wave 
energy  spectrum.  Then 


eu 

n(t)  = / cos  (uit-e(w))  /2S(u))dw 

0 


(C.l) 


= Re 


00 

j ei[0)t-e(u))]  /2S(w)du> 


(C.2) 


where  e(u>)  is  an  uniformly  distributed  random  variable  over  [0,2tt).  The 
expected  value  of  n(t)  is  then 

00 

En(t)  = E Re  J e*^  /2S(u)dw 


= Re 


00 

I 


eiWt  E{e"ie(u>)}  /2S(w)d^ 


The  characteristic  function  of  e(uO  is 


(t)  , ie(u)t,  = e27Tlt_^l 

Ve(u)  ~ Ei  ' 2irit 


(C.  3) 


which  implies 
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cp  (-t)  = E{e-ie(a))t} 
V £(w)  1 


. -2iTit 

1 - e 

2TTit 


Thus, 


(-1) 

e(o)) 


0 


E n(t)  = Re 


0 • /2S(w)dw  = 0 


(C.4) 


and  the  process  has  the  mean  value  zero. 

The  zero  mean  value  implies  that  the  covariance  function  of  n(t)  is 


cov  <n(t+x) ,n(t))  = e n(t+x)  n(t) 


= 1/2  Re 


+ 1/2  Re 


Now  for  w^0,  e(co), 
the  random  phase. 


00 

f ei[to(t+x)-0*t]  . E{e-i[e(ai)+e(0)]}  /2s(a,)2S(0)dwd0 
0 

00 

J ei[a)(t+T)+0*t]  . E|e-i[e(w)-e(0)]}  /2S(w)2s(0)du)d0 
0 (C.5) 

and  e(0)  are  statistically  independent,  by  definition  of 
Thus  the  characteristic  function  of  their  sum  is 


(C . 6) 


44 


E{e-iIe«o)+E<e)]}  , v(-«+e(e)  - 0,  we 


(C.  7) 


Similarly, 


'e(w)-e(6) 


«P(t)  • <p(t> 

ve(w)  v-e(9) 


COS  (2TTt)  - 1 

.2 


(C . 8) 


E{e-ileto)-e(e)],  , w(-«_£(e)  . 0,  we 


(C  .9) 


Thus,  applying  Equations  (C.7^  and  (C.9)  to  Equation  (C.5)  yields 


oo 

(n(t+T)n(t) ) = 1/2  Re  J*  E{e  2ie^w^}  2S(u))dw 


+ 1/2  Re  e1UT  E{l}2S(u))du) 


E{e-2i£(w)}  = *£»  = 0 


so  that 


COV  (n(t+T)n(t))  = 1/2  Re  j e 2S(ui)dai 

0 


c 

I 


iurr 


Oi 

'I 


cos  (u)T)S(w)du) 


(C.10) 


The  covariance  function  is  thus  a function  of  T and  is  independent  of  t. 
Hence,  the  wave  process  is  wide  sei.je  stationary.  We  write 


R(T)  = COV  (n(t+T)n(t))  = J S(w)  cos  (un:)dti) 

0 


0 

■J 


(C.ll) 


The  variance  of  the  process  is  then,  by  definition, 


a = R(0)  = j S(o))do) 
0 


f 


(C.12) 


The  probability  distribution  function  of  ri(t)  can  be  inferred  from 
the  definition  of  the  integral  of  Equation  (C.2)  by  expressing  the  integral 
as  a limiting  sum  and  properly  applying  the  central  limit  theorem. 

We  write 


n(t) 


Jt.  i.m. 

N 

Aw,  -*•  0 

k 


z cos  (w*t-ek)  |/2S((jj*)A(jJk 

k=l 


& • i • m • 

N 00 

Ao).  0 
k 


ii 


k=l 


(C. 13) 


46 


= £. i.m.  / \ 

N -*•  00 


(C  ’4) 


\ =/2S(u,k)Aa,k 


*k  " cos  (V"ek) 


(C. 15a) 


(C. 15b) 


Ak\ 


(C. 16) 


We  first  note  that  the  usual  central  limit  theorem  for  identically 

distributed  random  variables  is  not  applicable  to  the  sequence  (Yk). 

Fortunately,  a second  version  of  the  central  limit  theorem  can  be  used  to 

obtain  the  desired  result  provided  certain  regularity  conditions  are 

satisfied.  This  theorem  states  that  the  sum  of  a sequence  of  independently 

but  not  necessarily  identically  distributed  second  order  random  variables 

2 

with  means  pk  and  variances  ok  is  asymtotically  normally  distributed.  More 

exactly,  the  sum  converges  in  distribution  to  a normal  distribution.  A 

sufficient  condition  for  convergence  to  the  normal  law  in  this  case  is 
12 

given  by  a theorem  due  to  Liapounoff  which  requires  that  the  sequence  of 
third  absolute  moments  exist  for  the  random  variables  and  that 


£.i.m.  S = 0 
n 

n -*■  oo 


where 


F 


(C.17) 


It  is  not  obvious  by  inspection  that  the  condition  of  Equation  (C.17) 
holds  for  any  particular  application  of  the  central  limit  theorem.  The 
condition,  therefore,  must  be  verified  for  the  random  phase  model 
components . 

3 3 

Rice,  in  his  classic  paper,  proved  that  S -►  0 for  any  finite  range 

of  integration  in  the  random  phase  model  of  Equation  (C.14).  His  proof  can 

also  be  used  to  show  that  Sn  0 for  any  finite  frequency  range  (0,wc) . 

From  Equation  (C.14),  the  sequency  of  independent  random  variables  is 

defined  by 


Yk  = 


COS 


((\t~£k)’  k 1’2, 


(C. 18) 


The  variances  are  then 


°i=<4 


and  the  third  absolute  moments  are 


P I = E|Y.  |3 


E | cos3  (w*t-ek) | 3 


— a3 
3tt  \ 


(C.19) 


(C.20) 
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Then  from  Equations  (0.17),  (C.L9),  and  (0.20), 


Assume  Aui  •>  0,  u ► «>,  ? nAu>  » u>  . Then,  for  Aui  sufficient  Lv  small, 

r' 

1/3  ([' 

• I I S(w)dw 

' 0 

Thus,  since  nAu>  = u»  , we  have  S ► 0 as  n -►  m as  lone  as  at  is  finite. 

c n c 

However,  because  w can  be  chosen  as  large  as  desired,  we  say  that  the 
random  phase  model  from  EquaLion  (0.1)  converges  In  distribution  to  a 
normal  (Gaussian)  distribution  function. 

In  practice,  it  is  convenient  to  employ  probability  density  functions 
in  addition  to  probability  distribution  functions.  It  is  thus  of  interest 
to  know  whether  or  not  the  sequence  of  probability  density  functions  for 
the  random  phase  model  components  converges  to  a normal  density  function, 
von  Mlses  ^ has  given  the  following  sufficient  conditions  for  the 
component  random  variables: 
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1.  The  variances  satisfy  0 < < C < °“  Vk 

2.  The  third  central  moments  of  absolute  values  are  uniformly 
bounded. 

3.  The  character 1st ic  functions  <p^(t)  never  approach  the  value  1 at 
points  other  than  t = 0. 

A.  For  large  |t|,  the  <p^(t)  tend  to  zero  like  some  negative  power 
of  |t|,  and  uniformly  for  all  k. 

From  Equation  (C.18)  and  the  fact  that  the  £,  are  uniformly  distrib- 
uted, it  follows  that  probability  density  functions  for  the  are  given  by 


fkU) 


- A, 


x < A, 


k 


Thus 


o 


2 

k 


\ < 

2 S^cip  S(w)  < 00  for  every 


k 


and 


P 


3 

k 


4_  J 
3 it  k - 3n 


3/2 

S up  (S(w) ) 

U) 


for  every  k 


so  that  the  first  two  conditions  are  fulfilled.  Condition  3 is  also 
satisfied  because 


<Py 

yk 


dx 


<v> 


W 


where  J is  the  zeroth  order  Bessel  function  of  the  first  kind  which  equals 

° -1/2 
1 only  at  t “ 0.  For  large  |t|  J^A^t)  tends  to  zero  like  (tA^)  but 

not  necessarily  uniformly  in  k.  The  convergence  will  be  uniform  in  k if  it 

is  assumed  that  the  A^'s  are  bounded  below  by  a positive  constant.  The 

assumption  is  valid  for  any  frequency  range  [a,b]  satisfying  0 < a < b < 

Then,  since  a and  b can  be  selected  arbitrarily  close  to  zero  and  infinity, 

respectively,  we  say  that  condition  4 is  satisfied  for  the  random  phase 

model  and  conclude  that  the  sequence  of  density  functions  converges  to  a 

normal  density  for  each  t. 

14 

It  follows  by  a theorem  given  by  Lukacs  and  Laha  that  for  any  times 

t,  ...  t , the  distribution  of  (n(t.)  ...  r|(t  ))  = rf(t)  is  multivariate 
In  in 

normal.  Since  each  r)(t)  has  mean  zero  and  covariance  function  R(t),  rf(t) 
will  have  mean  vector  0 and  covariance  matrix  M = (R(t^-tj)).  The  finite 
dimensional  distributions  of  the  process  i"|(t)  depend  only  on  time 
differences  and  are  independent  of  t.  Therefore,  the  process  is  really 
strict  sense  stationary.  A strict  sense  Gaussian  process  is  ergodic  if 
and  only  if  its  spectral  distribution  function  has  no  points  of  dis- 
continuity as  expressed  by  Doob.^  This  condition  will  always  be  satisfied 
in  the  random  phase  model  for  the  following  reason.  We  have  assumed  that 
the  spectral  density  function  f (to)  exists  for  the  random  phase  model. 

Then,  from  Equation  (B.ll),  the  spectral  distribution  function  G^(io)  is 
given  by 


Gn«o) 


f n (9)d0 


S(0)d0 


Thus,  the  spectral  distribution  is  an  indefinite  integral.  Equivalently, 
this  means  that  G^(to)  is  an  absolutely  continuous  function,  and  thus,  it 
is  certainly  a continuous  function. 
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APPENDIX  D 

DIGITAL  GENERATION  OF  RANDOM  NUMBERS 

Random  numbers  are  required  for  the  digital  simulation  of  stochastic 
processes.  It  is  usually  sufficient  to  obtain  uniformly  distributed  random 
numbers  since  most  other  distributions  of  interest  can  be  obtained  as 
elementary  transformations  of  the  uniform  distribution. 

Ideally,  a random  number  represents  an  observation  from  some 
phenomenon  whose  outcomes  indicate  no  deterministic  regularity.  That  is, 
distinct  outcomes  may  be  different  even  when  observed  under  the  same 
conditions.  Thus  a sequence  of  random  numbers  will  never  repeat  itself 
and  it  can  be  identified  only  through  its  statistical  properties. 

In  the  past,  random  numbers  were  obtained  from  a number  of  diverse 
sources  such  as  roulette  wheels,  census  tables,  and  electronic  noise 
generators.  With  the  advent  of  high  speed  digital  computers  most  of  the 
random  number  sources  were  rejected  in  favor  of  simple  numerical 
algorithms.  Such  numerical  algorithms  produce  sequences  of  numbers  that 
satisfy  the  statistical  properties  of  randomness.  However,  the  same 
sequence  of  numbers  always  results  when  an  algorithm  is  exercised  twice 
under  exact  conditions.  Because  of  this  reproducibility  property, 
digitally  generated  numbers  are  called  "pseudorandom  numbers." 

There  are  various  algorithms  for  digitally  generating  pseudorandom 
numbers.  The  midsquare  method  appeared  quite  popular  at  one  time.  Its 
principle  is  as  follows.  Let  Xq  be  a p digit  number.  Its  square  gives 
roughly  2p  digits;  the  middle  p digits  are  taken  as  X^.  The  value  of  X^ 

is  then  squared  and  the  process  is  repeated.  Thus,  if  Xn  = 0.76, 

2 U 
Xq  = 0.5776  and  the  sequence  0.76,  0.77,  0.92,  0.46,  0.11,  0.12  ...  is 

obtained.  The  sequence  eventually  repeats  itself.  Also  for  certain 

choices  of  Xq,  such  sequences  converge  quickly  to  zero,  and  are  thus  biased 

toward  lower  values. 

Today,  most  computer  systems  possess  algorithms  for  generating 
pseudorandom  numbers  that  are  based  on  the  number  theory  concept  of 
binomial  congruence,  as  first  proposed  by  Lehmer.^^  Such  algorithms  are 
defined  by 
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PRECEDING  Pifls  hi  *Ky 


(D.l) 


X = KX  + b mod(2P) 
n+1  n 

where  p is  the  word  size  in  bits  for  the  employed  binary  machine.  The 
indicated  congruence  relation  means  that  if  the  difference  (Xn+^-KXn)  were 
divided  by  2P  the  remainder  would  be  b.  When  b = 0,  the  generator  is 
called  multiplicative;  otherwise  it  is  called  mixed.  Sequences  generated 
by  such  congruence  algorithms  also  have  finite  periods  (cycle  lengths) . 

The  constant  K is  chosen  to  maximize  the  cycle  length.  The  initial 
number,  X^,  is  usually  taken  as  any  odd  number. 

Random  numbers  generated  by  the  system  routine  on  CDC  6000  series 
computers  are  multiplicative  and  are  generated  according  to 

Xn+1  = (217+55205)  Xn  mod(l)  (D.2) 

where  mod(l)  implies  that  the  numbers  in  Equation  (D.l)  have  been  divided 
by  2P  to  obtain  pseudorandom  numbers  that  assume  values  between  zero  and 
one.  The  CDC  6000  computer  series  uses  sixty  bit  words  such  that  p = 60. 

As  implemented  by  CDC,  the  generator  is  designated  as  RANF  and  has  a full 
period  of  2^^. 

Given  a finite  sequence  of  numbers  generated  by  a numerical  algorithm, 
it  is  desirable  to  determine  whether  the  numbers  are  distinguishable  from 
a random  sample  or  from  a real  random  process.  The  standard  technique  for 
this  determination  has  been  the  use  of  standard  "goodness-of-f it"  tests. 
Canavos  has  performed  a comprehensive  set  of  goodness- to-fit  tests  in- 
volving the  CDC  6000  series  generator,  RANF.  His  tests  included  analyses 
of  runs,  serial  correlations,  moments,  chi-squares,  and  Kolmogorov- 
Smirnov  statistics.  Each  method  was  tested  against  a set  of  ten  sequences 
of  random  numbers  and  for  sample  sizes  of  200,  1,000,  and  10,000.  Al- 
though his  results  indicated  slightly  different  results  for  the  different 
tests,  the  overall  results  indicated  that  RANF  does  generate  uniformly 
distributed  pseudorandom  numbers  over  the  interval  (0,1).  To  further 
investigate  the  quality  of  numbers  obtained  from  RANF,  twelve  sequences 
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of  10,000  numbers  each  were  generated  on  the  DTNSRDC  CDC  6700  system  and 

chi-square  tests  were  performed  at  the  0.10  significance  level.  The  number 

of  class  intervals  k was  selected  optimally  in  accordance  with  the  procedure 
8 9 

of  Mann  and  Wald  as  detailed  by  Williams.  Thus,  165  intervals  were  em- 

ployed which  gave  164  degrees  of  freedom  for  the  tests  and  yielded  a 
2 

critical  value  of  y,,;  n in  = 187.56. 

^164, 0. 10 

Based  on  the  assumed  statistical  hypothesis,  no  more  than  one  out  of 
ten  samples  would  be  expected  to  yield  a chi-square  value  that  exceeded 
the  critical  value  if  the  samples  were  representative  of  an  uniform 
distribution.  The  chi-square  values  computed  from  the  12  samples  are 
shown  in  Table  6.  As  indicated  in  this  table,  none  of  the  computed  values 
exceeded  the  critical  value.  Thus,  it  was  concluded,  that  pseudorandom 
numbers  generated  using  RANF  exhibit  the  characteristics  of  an  uniformly 
distributed  random  variable. 
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PROGRAM  HAVESIM  ( INPUT  .OUTPUT  , T APE6*0UTPUT  , TAPE7) 

•••PROGRAM  TO  GENERATE  RANDOM  NAVE  PROCESS  USING  P IERSON-MOSKOHIT Z 
AMPLITUDE  SPECTRUM  BY  RAE  HURNITZ  CCOE  1524  MARCH  7,1973 
N IS  THE  NUMBER  OF  CASES 
CELT  IS  THE  TIME  INCREMENT  IN  SECONDS 
OELN  IS  THE  FREQUENCY  INCREMENT  IN  RACIANS/SECONO 
M IS  THE  NUMBER  OF  FREQUENCY  SUBDIVISIONS 
PM  IS  THE  TOTAL  NUMBER  OF  FREQUENCIES,  MM  = M ♦ 1 
NO  IS  THE  TOTAL  NUMBER  OF  POINTS  COMPUTED  FOR  HAVE  PROCESS 
CMEGMIN  IS  THE  MINIMUM  FREQUENCY  IN  RAOI ANS/SECOND 
CMEGMAX  IS  THE  MAXIMUM  FREQUENCY  IN  RADIANS/SECOND 
TINMIN  IS  THE  LENGTH  OF  THE  HAVE  PROCESS  IN  MINUTES 
TINSEC  IS  THE  LENGTH  OF  THE  HAVE  PROCESS  IN  SECONDS 
IFPLOT  = 1 PL CT  OUTPUT 
IFFLCT  = 2 DO  NOT  PLOT  OUTPUT 

COMMON  H (105), NST AR (1  OS) ,ST ARRAN < 105 ) ,S ( 1 05) , PMS ( 105) , 

♦ RANDOM (105)  ,XX(105)  ,HH( 105) , TIME (1805) ,H A VE ( 1605 ) , A, ARE APMS , 

♦ DELT ,DELH,P,MP,N,NN,NC, OMEGMAX,OMEGMIN,PI ,RMEAN, 

♦ TINMIN, TINSEC ,VARANCE,VK 

500  READ  100, N,M,MP, CELT, OME GM IN, OMEGMAX  , TINMIN 
READ  120, IFPLOT 

DELH  = (OMEGMAX  - OMEGMI N) /FLOAT ( M) 

TINSEC  = 60.  • TINMIN 
CNC  = TINSEC/OELT 
ONO  = CNO  ♦ 1. 

NO  = CNO 
PRINT  116 

PRINT  10 2, N, NO, MM, CMEGMIN, OME GM AX, OELH, DELT, TINMIN, TINSEC 
CALL  SPECTRA 
PRINT  116 
PRINT  104 

PRINT  106,  (H(I),HSTAR(I)  ,PMS(I) ,S (I) ,1*1,  MM) 

CALL  HAVES 
PRINT  116 

PRINT  108,(XX(I) ,1=1, MM) 

CALL  ARAMEAN 
DO  510  1*1, NO 
HAVE ( I ) = HAVE  (I)  - RMEAN 
510  CONTINUE 
PRINT  116 
PRINT  110 

PRINT  112,  (TIME (J) ,HAVE(J)  ,TIME(J*1),HAVE(J*1),TIME(J*2) , 

♦ HAVE (J+2) ,T IME (J*3) .HAVE ( J*3) , J*l,NO,4) 

CALL  ARAMEAN 

CALL  VAR 
PRINT  116 

PRINT  114,RMEAM,VARANCE«  ARE  ARMS 
HMIN  * 0.0 
HMAX  = 0.0 
00  515  1=1, NO 

IF (HA  VE ( I ) .GT*  HMAX)  HMAX  * HAVE(I) 

IF  ( HAVE ( I)  ,LT.  HMIN)  HMIN  * HAVE(I) 
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515  CONTINUE 
PRINT  116 

PRINT  117,NMIN,NMAX 
IF  flFPLOT.EQ. 1)535,540 
535  CALL  PLOTN 
540  IF  (N.EQ.l)  550,500 
100  FORMAT  1315, 4F10. 2) 

102  FORMAT  <20X, ’NUMBER  OF  INPUT  CASES  •» 

♦I3//20X, ’TOTAL  NUMBER  OF  POINTS  COMPUTED  FOR  NAVE  PROCESS  ’,15// 

♦ 20X, ’NUMBER  OF  FREQUENCIES  USEO  ’,1 7// 

♦ BOX, ’MINIMUM  FREQUENCY  IN  RAO I A NS /SECOND  *,F8.3// 

♦ 20 X , ’M  AX IMUP  FREQUENCY  IN  RAO I A NS /SECOND  *,F8.3// 

♦ 20 X,* FREQUENCY  INCREMENT  IN  RADIANS/SECOND  *,F6.3// 

♦20X,*TIME  INCREMENT  IN  SECONDS  *,F6.3// 

♦ 20X, ’LENGTH  OF  NAVE  PROCESS  IN  MINUTES  ’,F6.3// 

♦20X,*LENGTH  OF  NAVE  PROCESS  IN  SECONOS  *,F8.3//) 

104  FORMAT  < 2X , ’FREQUENCY  SUBDIVISION’, 5X, ’ACTUAL  FREQUENCY’, 5X, 

♦ ’SPECTRAL  OROINATE  0ASED’/48X, 

♦ ’ON  ACTUAL  FREQUENCY*//) 

106  FORMAT  C 8X, F6. 3,1 8X , F6. 3 , 16X ,F 8. 3 ,20X ,F8 . 3 ) 

110  FORMAT  (10X,*TINE’,7X,*WA VE*,1 OX, ’TIME*, 7X,’HAVE’, 10X, ’TIME’, 

♦ 7X,’NAVE*,10X,’TIME*,7X,  ’NAVE*//) 

112  FORMAT  <8X,F8.2,3X,F8.2,6X,F8.2,3X,F8.2,6X,F8.2,3X,F8.2,6X,Fe.2, 

♦ 3X  ,F  8 . 2) 

114  FORMAT  ( 20X, ’MEAN  OF  THE  GENERATED  NAVE  PROCESS 

♦ ’.F10.5//20X, ’VARIANCE  CF  THE  CCMPUTEO  NAVE  PROCESS 

♦ ’.F10.5//20X, 

♦ //20X»*AREA  BASEO  ON  THE  CONTINUOUS  FORM  CF  THE  PIERSON- 
♦MCSKOhlTZ  SPECTRUM  *,F10.5//) 

116  FORMAT  I 1H1 ) 

117  FORMAT  (2X,’HMIN  = ’ , F1C .5 , 10X ,’NM AX  = *,F10.5) 

120  FORMAT  115) 

550  STOP 
END 
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SUBROUTINE  SPECTRA 

THIS  SUBROUTINE  CALCULATES  THE  SPECTRAL  DENSITY  FUNCTION  USING 
THE  PIERSCN-HCSKONITZ  FORPULA 

RANFCCLM)  IS  A COC  EXTERNAL  FUNCTION  MHICH  RETURNS  VALUES  UNIFORHLY 
DISTRIBUTED  OVER  THE  RANGE <0  .1 ) .RANOOP  NUPBER  GENERATOR 
RANSET <A ) IS  A CDC  SUBROUTINE  MHICH  INITIALIZES  THE  VALUE ( A)  OF  THE 
RANDOM  NUMBER  GENERATOR 
A IS  A FLOATING  POINT  NUMBER 
PNS  IS  THE  PIERSON-MOSKOMITZ  SPECTRAL  ORDINATE 
BASED  ON  FREQUENCY  SUBDIVISION 
S IS  THE  PIERSCN-MOSKOHITZ  SPECTRAL  OROINATE 
EASED  ON  ACTUAL  FRECUENCY 
SCJ)  ARE  THE  SPECTRAL  OROINATE  VALUES 
STARRAN  IS  THE  RANOOM  NUMBER  OBTAINEO  FROP  RANF 
VK  IS  THE  HIND  SPEED  IN  KNOTS  FOR  THE  PIERSON-MOSKOMITZ  FORMULA 
M ARE  THE  FREQUENCY  SUBDIVISIONS 
MSTAR  ARE  THE  ACTUAL  FREQUENCIES 

COMMON  M <105) .MSTAR (105) , ST ARRAN < 105) ,S <1 0 5) , PMS (105) , 

♦ RANOCPC105) .XXC105) ,MN< 105).TIPE  <1 805 ) , M AVE < 1805 ) , A, AREAPPS , 

♦ delt.celm.p.mh.n.nn.no, cmegmax.opegpin.pi.rmean, 

♦ TINMIN.TINSEC.VARANCE.VK 
READ  100. VK 

PRINT  110, VK 
IOC  FORMAT  (F  10.2) 

110  FORMAT  <20X,*MINO  SPEEC  IN  KNOTS  * 

♦ Ffc.3//) 

M ( 1 ) * OMEGMIN 
CO  120  1*2. MP 
120  MCI)  * K»T-1)  ♦ OELM 
MSTARI1)  s OMEGMIN 
CALL  RANSET (3.  12) 

DO  130  1*2. MM 

ST  ARRAN  II)  s RANF (CUM) 

MSTAR  <I)  * MCI)  ♦ STARRANCI)  * OELM 
130  CONTINUE 

A 3 .00810  * <32.174)**2 
B 3 .74  • <3  2.  174/CVKM.6878)  )**4 
PI  * 3.1415927 
CO  140  J=1 »MM 

SCJ)  = A/CMSTAR<J)**5)*EXPC-B/<MSTAR<J)**4)) 

PMSCJ)  3 A/CMCJ)**5)*EXP C-E/ CM C J) **4) ) 

140  CONTINUE 
RETURN 
END 
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SUBROUTINE  NAVES 

THIS  SUBROUTINE  COMPUTES  THE  RANDCH  NAVE  PROCESS 
RANOOM  IS  A RANDOM  NUMBER  BETNEEN  0 AND  1 

RAhF(OUM)  IS  A COC  EXTERNAL  FUNCTION  NHICH  RETURNS  VALUES  UNIFORMLY 
0 ISTR I6UTE0  OVER  THE  RANGE  10*1). RANOOM  NUMBER  GENERATOR 
RANSETU)  IS  A COC  SUBROUTINE  NHICH  INITIALIZES  THE  VALUE  (A)  OF  THE 
RANCOM  NUMBER  GENERATOR 
TIME  IS  TIME  IN  SECONOS 

NAVE  IS  THE  NAVE  PROCESS  AS  A FUNCTION  OF  TIME 
XX  IS  THE  RANOCM  PHASE  ANGLE 

COMMON  N (105) * NSTAR(105) *ST ARRAN (105) ,S(105) ,PMS (105) , 

♦ RANOCM (105) *XX (105) ,NN( 105), TIME (1805), NAVE (1005), A, ARE ARMS, 

4 oelt,oeln,m,mm,n,nn,no,omegmax,omeghin,pi,rmean, 

♦ TINMIN.TINSEC  .VARANCE  ,VK 
TIME  BEGINS  AT  ZERO 

TIME  (II  * 0.0 
00  20  JJ=2,NC 

20  T I ME  ( J J I = TIME(JJ-l)  ♦ OELT 
CALL  RANSET  (7.111 
00  25  1*1, MM 
RANOOM (I)  = RANF(OUM) 

25  XX(I>  = RANOCM  (I)  * 2.  * PI 
00  30  X* 1 ,NO 
30  NAVE (XI  * 0.0 
00  40  J=l, NO 
00  35  1=1, MM 

WH(I>  = SCRT (2.*S(II*0ELN) *COS ( NSTAR ( I )*T IME ( J ) - XX(I)) 

NAVE ( J)  = NAVE  (Jl  ♦ NN (I  I 
35  CONTINUE 
40  CONTINUE 
RETURN 
END 
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SUBROUTINE  ARAMEAN 

AREAPMS  IS  THE  AREA  BASEO  OH  THE  COHTIHUOUS  FORH  OF  THE 
PIERSON-MOSKOMITZ  SPECTRUM 

THIS  SUBROUTINE  CALCULATES  THE  MEAN  VALUE  OF  THE  GENERATED  HAVE  PROCESS 

RMEAN  IS  THE  T IME- AVER AGE  OF  THE  HAVE  PROCESS 

COHMOK  H(105),HSTAR(10S) ,ST ARRAN (105) ,S< 1 051 , PMS ( 1051 , 

♦ RANDOM! 105) , XX (105) , HH( 105) , TIME (1005) . HA VE ( 1805 ) . A, ARE APHS . 

♦ OELT.CELH.M.MN.N.t.N.NO.ONEGHAX.OMEGMIN.PI.RNEAN, 

♦ T INHIN.  T INSEC  .VARANCE  »VK 
COM  « . 00810*  (32. 174)**2 

COM2  s .74  * (32.174/(VK*1.6878))**4 

AREAPMS  = CCN1/I4.*C0N2) * (EXP (-C0N2/0MEGMAX**4)- 

♦ EXP(-C0N2/0MEGMIN**4) ) 

HMEAN  = 0.0 

FNO  = FLOAT(NO) 

00  60  1=1. NO 

60  HMEAN  * HAVE (I ) ♦ HMEAN 
RMEAN  = HMEAN/RNO 
RETURH 
END 


SUBROUTINE  VAR 

THIS  SUBROUTINE  CALCULATES  THE  VARIANCE  OF  THE  COHPUTEC 
NAVE  PROCESS 

CONNOR  H(105)»VSTAR(105) , ST ARRAN f 105) ,S<105) ,PMS<105) , 

♦ RANDOM 105) ,XX (1 05) , NN ( 105 ) , TIME (1 805) , MA VE 1 1805 ) , A. AREAPNS , 

♦ OELT,OELM,M,H)>,NtRN'NC'CMEGMAX,OMEGMINtPI,RHEAN, 

♦ TINMIN,TINSEC,VARANCE,VK 
VARIAR  a 0.0 

RNO  * FLOAT  I NO ) 

CO  70  1=1, NO 

70  VARIAR  = VARIAN  ♦ HAVE(I)  * MAVE(I) 

VARANCE  = VARI AN/RNO 

RETURN 

ENO 


SUBROUTINE  PLOTH 

C ONION  W (1 051, HST AR(105)  , STARR AN < 1051 ,S( 1051 ,PMS(105) , 

♦ RANDOM (105) , X X ( 105) ,HH(105) ,T IME t 1805 ), NAVE C 1805) , A, ARE  ARMS, 

♦ DEL  T,  HELM  ,N,NH , N , NN, NO , OMEGMA X, CMEGNTN . PI , RME AN , 

♦ tinnin,tinsec .varance »vk 

C PLCT  PCUTINE  FOR  CALCCMP  936 

CALL  PLOTS(0.0«0.0.7) 

CALL  PLOT  (0.0,0. 0,-3) 

CALL  SCALEUINE, 18.0, 1801,1) 

CALL  SCALE(MAVE, 4. 0,1801,1) 

CALL  AXIS(C.O,0. C,15HTIME  IN  SECCNCS, -15,18. 0,0.0, TIME  (1802)  , 

♦ TIME (1803)) 

CALL  AXIS  (0.0  ,0.0,14HHAVE  ELEVATION, 14, 4. 0,90.0, HAVE  ( 1802)  , 

♦ WAVE ( 1803) > 

CALL  LINE  (TIME, HAVE, NO, 1,0,0) 

CALL  SVMPOL(2. 0,8.0,. 20, 12HTIME  HIS  TORY , 0 , 1 2) 

CALL  PLOT (0.0,0. U, -3) 

CALL  PLOT(0.0,0. (,999) 

PRINT  650 

650  FORMAT  (20H1FLCT  IS  FINISHED  ) 

RETURN 

ENO 
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